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We test three different approaches, based on quantum Monte Carlo simulations, for computing 
the velocity c of triplet excitations in antiferromagnets. We consider the standard S = 1/2 one- 
and two-dimensional Heisenberg models, as well as a bilayer Heisenberg model at its critical point. 

Computing correlation functions in imaginary time and using their long-time behavior, we extract 
the lowest excitation energy versus momentum using improved fitting procedures and a generalized 
moment method. The velocity is then obtained from the dispersion relation. We also exploit 
winding numbers to define a cubic space-time geometry, where the velocity is obtained as the 
ratio of the spatial and temporal lengths of the system when all winding number fluctuations are 
equal. The two methods give consistent results for both ordered and critical systems, but the 
winding-number estimator is more precise. For the Heisenberg chain, we accurately reproduce 
the exactly known velocity. For the two-dimensional Heisenberg model, our results are consistent 
with other recent calculations, but with an improved statistical precision; c = 1.65847(4). We 
also use the hydrodynamic relation c? = ps/X-c{<l 0) between c, the spin stiffness ps, and the 
transversal susceptibility xxi using the smallest non-zero momentum q = 2'itIL. This method 
also is well controlled in two dimensions, but the cubic criterion for winding numbers delivers 
better numerical precision. In one dimension the hydrodynamic relation is affected by logarithmic 
corrections which make accurate extrapolations difficult. As an application of the winding-number 
method, for the quantum-critical bilayer model our high-precision determination of the velocity 
enables us to quantitatively test, at an unprecedented level, field-theoretic predictions for low- 
temperature scaling forms where c enters. We find agreement to within 3% with 1/A expansions 
for the coefficients of the leading susceptibility and specific heat forms; x ^ and C ~ hT^. 

PACS numbers: 75.10.Jm, 75.40.Mg, 75.30.Ds, 75.40.Cx 


I. INTRODUCTION 

Ordered quantum antiferromagnets exhibit linear dis¬ 
persions of their elementary spin-wave (magnon) exci¬ 
tations and the associated velocity c is an important 
parameter characterizing such systems. The prototyp¬ 
ical two-dimensional Heisenberg model has an ordered 
ground state and its low-energy magnon spectrum is well 
described by spin wave theoryji When 1/S corrections 
are properly taken into account, the velocity and other 
properties computed within this approximation for the 
most extreme (and interesting) case of the spin S =1/2 
model agrees to within « 1% with results of quantum 
Monte Carlo (QMC) calculations, which can be consid¬ 
ered exact to within statistical errors if sufficiently large 
lattices are used for extrapolations to the thermodynamic 
limiti^^— This good agreement can be traced to the fact 
that the ground state is strongly ordered, the sublattice 
magnetization being reduced by quantum fluctuations by 
only about 40% from the classical value. Upon introduc¬ 
ing other interactions which enhance the quantum fluctu¬ 
ations and suppress the order, the quantitative predictive 
power of spin-wave calculations rapidly deteriorates and 
the quantum fluctuations have to be treated in more so¬ 
phisticated waysi^^— An extreme case is when a system 
is driven to criticality. The low-energy critical excita¬ 
tions are still linearly dispersing but the corresponding 


velocity cannot be reliably calculated in any simple the¬ 
oretical manner. Lastly, quantum-disordered antiferro¬ 
magnets also have propagating triplet excitations, which 
are gapped and often called triplons. 

In this paper we will discuss three very different ways 
to extract the velocity of the elementary excitations of 
quantum spin models based on ground-state projector 
and finite-temperature QMC simulations. We consider 
one-dimensional (ID) and two-dimensional (2D) ordered, 
disordered and critical quantum antiferromagnets. Using 
imaginary-time dependent spin correlation functions, the 
long-time behavior computed at different momenta con¬ 
tain information on the dispersion relation, from which 
the velocity can be extracted if the limits of the system 
size going to infinity and the momentum going to zero are 
treated correctly. One can also in some cases extract the 
velocity in a simpler, indirect way using winding num¬ 
ber fluctuations in the space-time representation of the 
system sampled in the QMC calculations We will 
develop stable procedure based on these two approaches 
and compare the results for a number of different cases. 
In addition, we also test the well known hydrodynamic 
relationship = Ps/X-L for antiferromagnetic state/^ 
where X-L is the transversal magnetic susceptibility and 
Ps the spin stiffness. 

We first successfully test the methods on both ID 
and 2D systems for which the velocity (of spinons and 
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spin waves, respectively) is previously known, and there¬ 
after study a bilayer system, both in its paramagnetic 
phase and at its quantum-critical point. In the lat¬ 
ter case we subsequently use our high-precision estimate 
of c to investigate in detail, to our knowledge at an 
unprecedented level of control, the reliability of finite- 
temperature quantum-critical scaling forms for the mag¬ 
netic susceptibility and the specific heatj ^^d^ 

Before turning to the calculations and results, in 
Sec. ITU we first provide some more background on the 
utility of carrying out precise determinations of the ve¬ 
locity in quantum spin systems, focusing on quantum- 
criticality in dimerized antiferromagnets. In II Bl we pro¬ 
vide a brief summary of the different calculations to be 
presented and in EC] we outline the organization of the 
rest of the paper. 


A. Effective low-energy descriptions 

Quantum field theories are often used to describe uni¬ 
versal low-energy properties of quantum magnetsi ^^d^ 
Numerical techniques, such as QMC simulations, can be 
used to extract system-dependent parameters appearing 
in various predicted forms of physical quantities at low 
energy. Such an approach of combining quantum field 
theory and numerics has been established over the past 
several years for certain types of quantum phase tran¬ 
sitions in 2D systemsj ^'^d^ Most well studied are tran¬ 
sitions in dimerized models, where for weak inter-dimer 
coupling there is a tendency to singlet formation on the 
dimers, leading to a quantum phase transition into a 
quantum paramagnet at a critical ratio of the inter- 
and intra-dimer couplings ° d A similar transition 
(of the same universality class) also take place if the 
dimers are replaced by some other unit cell of an even 
number of spins on which the single-cell ground state is 
a singletj ^^d^ The many studies of these systems have 
shown rather convincingly that the phase transition is 
in the 3D 0(3) universality class, in agreement with 
field-theory descriptions based on the non-linear sigma 
model] Other properties associated with quantum- 
criticality are also well captured by the field theory^ 
e.g., the uniform magnetic susceptibility x is linear in 
temperature at the critical coupling ratio, x = aT/c^, 
and the specific heat grows quadratically; C = hT'^ je? at 
low T. Perhaps surprisingly, however, the values of the 
prefactors a and b have still not been tested quantita¬ 
tively in a completely unbiased manner)^ This may be 
largely because, as indicated above, they depend on the 
velocity c of the critically damped spin waves (and on no 
other low-energy parameter), but this parameter has not 
yet been independently calculated to high precision for 
model systems (by first-principles methods not depend¬ 
ing on other field-theory predictions). As a demonstra¬ 
tion of the value of determining c to high precision, in 
this paper we will also provide a test case of a detailed 
comparison with field theory predictions for one of the 


prototypical dimerized models; the Heisenberg bilayer. 


B. Technical and Physics Objectives 

The first aim of the present paper is to systematically 
test three completely different ways of extracting the ve¬ 
locity of elementary excitations based on QMC calcula¬ 
tions in the following ways: (i) Using the momentum de¬ 
pendent gaps extracted from imaginary-time dependent 
spin correlation functions, (ii) By monitoring spatial and 
temporal winding number fluctuations, which are propor¬ 
tional to the spin stiffness and the uniform susceptibility, 
respectively, and adjusting the space-time aspect ratio 
L/fi such that these fluctuations are equal. At this spe¬ 
cial inverse temperature fd*{L) the ratio of the spatial and 
temporal lengths L/(3* should equal c.— (in) Exploit¬ 
ing the hydrodynamic relation = Psjx^^ where one 
has pay attention to the fact that xi. 0 when the tem¬ 
perature r —^ 0 in a finite system, while the spin stiffness 
remains nonzero. We discuss an approach to circumvent 
this problem based on X-l('?) where the momentum q is 
small but non-zero. 

The method (i) is in principle very direct, being con¬ 
nected to the fundamental definition of c in the dis¬ 
persion relation of the lowest-energy excitation versus 
momentum. However, the precise determination of the 
momentum-dependent gap is in practice complicated by 
the presence of a continuum of excitations above the low¬ 
est energy. In some cases, as we will show, one also has 
to take great care with the way the thermodynamic limit 
and zero-momentum limits are taken. The method (ii) 
is rather simple and the only uncertainty is introduced 
by a final extrapolation of the finite-size velocity defi¬ 
nition Ljfd* to the thermodynamic limit. However, as 
far as we are aware, the correctness of this approach has 
not been formally proven, except for the case of a long- 
range ordered antiferromagnet, and the method has not 
been widely usedj ^^d^ We here confirm that the method 
continues to give the correct velocity also when the sys¬ 
tem is critical, even in the case of the S' = 1/2 Heisen¬ 
berg chain, where the low-energy excitations are not even 
spin-waves but deconfined spin-1/2 carrying topological 
defects (“spinons”). The method (iii) based on general¬ 
ized hydrodynamics is simple to apply and we will argue 
that it as well continues to work also in critical systems. 

Comparing methods (i)-(iii) for the Heisenberg chain 
as well as for the well-studied 2D Heisenberg model gives 
us insights into how to best apply the methods in prac¬ 
tice. After these tests, we study the bilayer Heisenberg 
model at its quantum-critical point and in its paramag¬ 
netic regime, using methods (i) and (ii) Also in this case 
we find good agreement between the two methods at the 
critical point, but one has to be more careful when defin¬ 
ing the velocity based on gaps for finite systems, because 
of a slower convergence of the gaps to their infinite-size 
values than in the ordered state. We also explicitely show 
that the winding number estimator gives an incorrect ve- 
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locity in the gapped paramagnetic phase. 

The second aim of the paper connects to the low-energy 
filed-theory description discussed in Sec. II Al — to com¬ 
pute a high-precision value for c for the quantum-critical 
Heisenberg bilayer model, and to use this value to reli¬ 
ably test the field-theory predictions for x(T) and C{T). 
While many such tests have been performed in the past, 
on the bilayer— as well as other— quantum-critical 
2D spin systems, the past studies were not able to com¬ 
pletely quantitatively test the level of agreement with 
the existing large-field theory predictions, because an 
independent, unbiased determination of c was lacking. 
This obstacle is overcome by the reliable calculation of c 
in this paper. 


triplet gap (with respect to the ground state energy) as 
a function of momentum k. Applying the imaginary- 
time evolution operator to a trial state leads to 

the ground state when /3 —?► oo. Alternatively, one can 
use a high power iJ™ of the Hamiltonian and reach the 
ground state when m —> oo. In practice QMC simulation 
methods based on these two operators are very similar, 
but the exponential form allows more direct access to the 
standard imaginary time. As we will show, this connec¬ 
tion is not needed if only the excitation energies are of 
interest, and then one can use the slightly faster power 
method. We will use both approaches here, with two 
different ways of analyzing the correlation functions. 


C. Outline of the Paper 

The rest of the paper is organized as follows: In 
Sec.ini we discuss calculations of imaginary-time corre¬ 
lation functions within a ground-state projector QMC 
method formulated in the basis of valence bonds. We 
discuss our methods to extract the lowest momentum- 
dependent gaps from the long-time behavior of these cor¬ 
relations. In addition to fitting a leading exponential 
decay with additional corrections to account for higher 
excitations, we also present a generalization to ground- 
state projection QMC data of a systematic high-order 
moment approach recently introduced for use with finite- 
temperature QMC simulations— We compare our results 
with exact solutions for small Heisenberg chains as well 
as with the rigorously known velocity of this model in 
the thermodynamic limit. In Sec. IHII we discuss the de¬ 
termination of c based on winding numbers and demon¬ 
strate the method using the Heisenberg chain. These 
finite-temperature calculations are carried out with the 
stochastic series expansion (SSE) QMC method. In 
Sec.HV] we discuss details of the hydrodynamic relation¬ 
ship generalized to finite system size and test it using 
SSE calculations for the Heisenberg chain. In Sec. IVl we 
present further tests of all three methods for the stan¬ 
dard 2D Heisenberg model. We also show that the scal¬ 
ing of the triplet gap at momentum k = (tt, tt) is en¬ 
tirely consistent with quantum rotor-states carrying spin 
5 = 1. In Sec. ED we discuss results for the critical 
and disordered bilayer systems, including the compari¬ 
son of c determined using both methods. The analysis 
of the quantum-critical susceptibility and specific heat 
computed using large-scale SSE calculations is presented 
in Sec. lVHI We briefly summarize and discuss our results 
further in Sec. EUD 


II. MOMENTUM-DEPENDENT SPIN GAPS 

Here we discuss how to use zero temperature (T = 0) 
projector QMC methods to calculate imaginary-time cor¬ 
relation functions which are then used to extract the 


A. H-power projection 


This QMC approach is based on projection with a suf¬ 
ficiently high power of the hamiltonian on a trial state 
I'f/'t)- The process can be conveniently expressed in the 
energy eigenbasis of H, leading to the following expres¬ 
sion for the dependence on m: 


{-Hr\^t)=co{-E^y 
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Here |n), n = O,-- - ,A — 1 are the energy eigenstates 
of E[ and Eq is assumed to be negative, with its abso¬ 
lute value |Eo| being the largest in magnitude of all the 
energies (which can always be achieved by adding a suit¬ 
able constant to the hamiltonian). Then, if the expan¬ 
sion coefficient cq ^ 0 (which in practice is essentially 
guaranteed for any reasonable choisce of ipt)), we have 
oc |0) for large m, and the expectation value 
of an operator O at T = 0 can be written as 


( 0 ) = 




( 2 ) 


For the SU(2) invariant spin models considered in this 
paper, this form of the expectation value can be evalu¬ 
ated by importance-sampling, using a formulation of the 
projector QMC method in the non-orthogonal valence 
bond basis— An efficient way of sampling the contribu¬ 
tions to (—iJ)^"*, very similar to “operator-loop” updates 
developed within the finite-T SSE method^S can be for¬ 
mulated using loop updates in a combined space of spin 
components (5^) and valence bonds— The trial state is 
also expressed using valence bonds, in the form of an 
amplitude-product state— The details of the state (the 
form of the amplitudes) are not important, as the good 
convergence to the ground state is observed even if the 
state is not optimized— 

We will use this valence-bond variant of the projector 
QMC method for computing appropriate imaginary-time 
dependent correlation functions. For technical details on 
the sampling methods we refer to Ref. [^. Below we will 
focus on the definition of the correlation functions we 
study and how we process them to extract the velocity. 
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1. Imaginary-time correlations 


We consider correlation functions of the following form: 


CA{t) 


{0\A^-HYA\0) 

m-H)m 
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|(n|^|0)p, (3) 


where t is an integer which can be related to imaginary 
time^ and we will loosely refer to it by this term. More 
precisely, t/N^ where N is the system size, is propor¬ 
tional to imaginary time t in the sense of the standard 
Schrodinger time evolution operator e“'^^, as we will ex¬ 
plicitly show below. From C^(f), we further define 


QAit) = 


CA{t)-mA\0)Y 

C^(0)-|(0|^|0)P’ 


( 4 ) 


and note that Qa( 0) = 1 and QaY —>■ oo) —>■ 0. For large 
t, we have 


( i(i|A|o)p ^ 

yo|AU|o)-|(o|A|o)pj \ey) ’ ^ ^ 

where A = i?i — i?o is the energy gap between the first 
excited state (connected to the ground state by the op¬ 
erator A) and the ground state. To directly show the 
relationship between t and imaginary time, we can intro¬ 
duce r = t/lEol, then write Eq = Ncq, and in the limit 
of large N have 

/ A \ Af|eo|T 

Qa{t = t/\Eo\) (X ^ (6) 

which is the familiar form of the asymptotic decay of 
an imaginary-time dependent correlation function. Thus, 
we have shown that, indeed, r oc t/N. We will not ex¬ 
plicitly need to make use of the relationship between t 
and r here, however, and we will continue to use t as the 
“time” parameter with the iJ-power method. 

We can appropriately choose the operator A so that 
it excites the ground state |0) into a state with desired 
quantum numbers. The ground state of the unfrustrated 
hamiltonians considered here are total spin singlets with 
momentum k = 0 on a finite lattice with an even num¬ 
ber of spins and periodic boundary conditions (in one 
dimension only when the number N of sites is a multiple 
of four—for other even N the momentum is tt)^ Since 
we are interested in triplet excitations, we can use the 
simple Fourier-transformed spin operator, 

A(k)=5^(k) = ^5^e*>^--5^(r), (7) 


to create an S' = 1 state with momentum k and = 0 
when acting on the ground state. Thus, the following 
imaginary-time correlation function 


Cu(t) 


(0|7l(-k)(-g)*A(k)|0) 

(■iPt\{-H)^^-P-*Ai-k){-HyAi-H)P\Yt) 


allows us to directly measure the triplet excitation gap 
A(k) as a function of the momentum. The second line in 
the above equation explicitly shows the form used in the 
projector QMC calculations, where both 2m — p — t and 
p are assumed to be large enough to achieve projection 
to the ground state for the system sizes considered. 

In practice, we will use projection powers m = Idp^N 
and, to achieve good ground-state convergence, require 
that 2m — p—t and p both are larger than 15pj.A, where 
typically p^ = 8 (16 or higher in some cases). These 
choices are motivated by extensive tests indicating that 
no detectable systematical errors remain. The values of t 
are restricted to be multiples of A/4 in our simulations. 
Since the correlation functions at the different t values 
are measured in the same simulation, the C'k(t) data are 
correlated and measuring at shorter t intervals does not 
significantly increase the amount of statistical informa¬ 
tion in the data set. 


2. Extracting the gap 


We here use two different ways to extract the lowest 
triplet gap from the correlation function Ck(t). Since 
(0|A(k)|0) = 0 for k ^ 0, Eq. (g]) reduces to 


Qk(t) 


Ck(t) 
Ck(o) ■ 


(9) 


In principle, the gap can be extracted by monitoring the 
long-time behavior, given by the form ([S])- A systematic 
way to extract the gap without performing any curve fits 
is to consider the ratio of Qk{t) at two different times 
separated by some interval; e.g., by A/4 operations: 


i?k(i) = \Eo\ 1- 


Qk{t + A/4) 


1 4/ArN 


Quit) 


( 10 ) 


Note that Rk{t) A(k) when t is large enough, which 
follows from the long-time behavior of (5k(i) in Eq. ([5|). 
However, in our QMC calculations, it is not always pos¬ 
sible to reach perfect convergence of this gap estimate for 
all k, because the relative statistical errors often become 
too large already for moderately large t. This problem is 
related to the existence of a continuum (for large A) of 
states above the gap, due to which the pure exponential 
decay cannot be easily observed in practice. 

We thus use another method to estimate the value of 
the gap based on the entire available set of correlation 
functions. This scheme is more reliable than the ratio 
scheme when data (with small relative errors) are not 
available for large values of imaginary time. It is clear 
from Eq. m that one can define a positive-definite spec¬ 
tral function Hk(w) to fit the normalized imaginary-time 
correlation (5k(i) as 

Qk(i) = ^ da;Ak(w) (^1 - . 


( 11 ) 
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This is just the analogue for the i/-power evolution of 
the standard form, 

pOO 

Gk(r)= / dw5u(a;)e-“", (12) 

Jo 

relating the imaginary-time Schrddinger evolved correla¬ 
tion function 


G^{T) = {0\S\{T)Sm\0). ( 13 ) 

to the dynamic spin structure factor 5'k(a;). In Eq. dill) uj 
cannot exceed | Eq |, following from the fact that we have 
ensured that Eq is the eigenvalue with the largest magni¬ 
tude. In practice, as in the standard dynamic spin struc¬ 
ture factor in m, the actual dominant spectral weight 
will be concentrated only to within a window of order J. 

For any finite system, is a sum of (5-functions and 
this can be replaced by a continuum starting at w = A(k) 
for a large system (or, in some cases, there is an iso¬ 
lated ^-function at the gap, following by a second gap 
and then a continuum). With or Qk{t) computed 

using QMC calculations, the respective relations (ED or 
ED can in principle be inverted using numerical ana¬ 
lytic continuation. This procedure is very challenging, 
however, and it is not easy to extract the gap precisely 
with conventional methods such as the Maximum En¬ 
tropy method^ though one can extract the main domi¬ 
nant spectral features (and we note that progress in this 
regard has been made very recently^). Here our goal is 
merely to extract the gap, and instead of trying to repro¬ 
duce the full spectral function we model the excitations 
by just a small number of (5-functions. With the precision 
of typical QMC data, Qk(I) can be normally fitted very 
well with just a few (5-functions (typically 3 to 5) over 
the full range of accessible times t. With this procedure, 
we expect the location uji of the lowest gap to accurately 
reproduce A(k), while the higher (5-functions represent 
approximately the contributions of the continuum. In 
this fitting procedure, the extracted lo\ is to some extent 
affected by contributions of the higher states but does 
not change significantly when increasing the number of 
(5-functions. We can therefore quite reliably extract the 
lowest gap, but not higher ones unless they are separated 
by significant subsequent gaps (which is not expected in 
the cases of interest here, except well inside the quantum 
paramagnetic state of the bilayer model). 

Given a set of n ^-functions at energies uji with asso¬ 
ciated amplitudes Ai normalized so that Ai = 1, one 
can compute the associated time dependent correlation 
function in analogy with Eq. ED as 


Qk(t)= 



(14) 


Now denoting the corresponding QMC-computed func¬ 
tion by and their statistical error by (Tt, the good¬ 

ness of the fit is quantified in the standard way by 


based on a set of Nt time points {t}: 

^' = ^EA(Qk(t)-Qk(i))'- (15) 


We here use a uniform grid of time points with sepa¬ 
ration N/4, operations, t = 7V/4, iV/2,..., up to a point 
imax = lVt(5V/4) where the relative statistical error of 
Q\s,{t) exceeds 5%. The choice of cut-off is not very im¬ 
portant as, in any case, the noisy data at very large t will 
not affect the fit from the definition of . 

For the extracted gap to be reliable, the contribution 
of the lowest (5-function to the fit must be significant at 
the longest times included. To monitor this long-time 
weight, we compute the relative contribution of the low¬ 
est (5-function, denoted by Ai{t), at the time tmax in¬ 
cluded in the fit: 


Ai{t 


Hi(l-a;i/|£;o|)*-"- 


(16) 


This quantity should approach 1 for t —>■ cx) if the lowest 
(5-function is at the gap. It is close to 1 in all the fits 
reported here, indicating stable extraction of uji. 

The statistical error of the extracted gap is estimated 
using a bootstrap error analysis. With the underlying 
QMC data for the correlation function Gk(t) saved as M 
bin averages (with typically M ^ 100— 1000), bootstrap 
averages are constructed by selecting M bins at random 
(i.e., allowing for the same bin to be selected multiple 
times). The above fitting procedures are then carried 
out repeatedly for a large number of these samples, and 
the standard deviation of the estimates is the statistical 
error of the gap in our procedure. 

As already mentioned, the fluctuations of the QMC 
data at different times are significantly correlated since 
these are measured in the same simulation. A statis¬ 
tically correct treatment of the data would require the 
inclusion of the full covariance matrix (instead of just 
its diagonal elements) in the definition of However, 
much of the covariance is already removed when the time- 
correlations are normalized [by the denominator Eq. (|9D] , 
because the errors are correlated primarily by overall 
fluctuations in the normalization. Based on test cases, 
including ones reported below, to obtain fully reliable re¬ 
sults it is sufficient to use only the diagonal elements of 
the covariance matrix and define x^ as in ED- 


3. Tests on the Heisenberg chain 

We here illustrate the gap extraction method described 
in the previous subsection using the example of the 
S' = 1/2 Heisenberg spin chain with periodic boundary 
conditions, where spins interact with nearest neighbor 
exchange constant J = 1; 

L 

h = jJ2s.-s,. 

2=1 


( 17 ) 
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Relative error 


FIG. 1: (Color online) (a) Correlation ratio i? 7 r(t), Eq. (IIOII . 
shown for two Heisenberg chain sizes, L = 12 and 24. The 
lines are the exact values of the triplet gap at this wave-vector 
(with the numerical values indicated as well), (b) The nor¬ 
malized correlation function Q-,,{£) for L = 12 and L = 24. 
(c) The distribution of the relative error of the gap A(7r) ex¬ 
tracted by htting the data for L = 24 in (b) using three 5- 
functions. The histogram was generated from a large number 
of bootstrap samples of the QMC data. The relative error is 
calculated as (A — Ae)/Ae where Ae is the exact value of the 
gap. The exact value (corresponding to relative error 0) is 
seen to fall within a standard deviation of the distribution of 
the bootstrap samples. 

First, we compare the results of our numerics for the 
lowest triplet gap at fc = tt for chain sizes L = 12 and 
24 with exact diagonalization (Lanczos-method) results. 
As can be seen in Figllja), the quantity R-jrit) defined 
in Eq. (flUl) indeed converges to the correct gap value in 
both the cases. Fig [T](b) shows the normalized imagi- 



FIG. 2: (Color online) Convergence of the correlation ratio 
Rk(t) to the gap A(k) as a function of imaginary time t, 
shown at fc = TT and k = 27r/I/ for a Heisenberg chain of 
length L = 64. The inset shows Rk{t) at k = tt + 2tt/L for a 
longer chain, L — 200, where the gap extraction using a fit of 
Qk(t) to 5 5-functions works demonstrably better, giving the 
result indicated below the horizontal line, with the thickness 
of the line representing the average gap ± one error bar. 

nary time correlation function Q 7 r(t) for the two system 
sizes, and Fig[IJc) shows the distribution of the gap error 
obtained using a bootstrap method (i.e., the simulation 
data stored as M bin averages are re-sampled by select¬ 
ing M bins at random, and the fitting procedures are 
carried out for each such averaged data set) for L = 24, 
using a fit to three ^-functions. 

This analysis show that the gap obtained from the fit 
agrees statistically with the exact result, which here is 
within a standard deviation of the distribution obtained 
in the bootstrapping procedure, and the distribution it¬ 
self closely matches a normal distribution. Thus, even 
with only three 5-functions in the spectrum (which is 
clearly a much smaller number than what is contained 
in the full spectrum) we detect no systematic errors in¬ 
troduced by the fitted functional form, supporting our 
assertion that the simplified description of the spectrum 
does not significantly affect the location of the lowest 5- 
function (the gap). The number n of 5-functions that 
should be included in a given case depends on the statis¬ 
tical errors of the QMC data and the actual form of the 
spectral continuum. To determine n, we monitor as a 
function of increasing n and stop when no improvements 
in the fit are observed. 

For the Heisenberg chain, the lower edge of the spec¬ 
trum of triplet excitations is known rigorously based on 
the Bethe Ansatz solutioni^ For momentum /c —>■ 0 and 
/c —> TT the spectrum is linear with velocity c = 7r/2. 
For an infinite chain the triplets are degenerate with sin¬ 
glet excitations, due to the fact that the excitations con¬ 
sist of essentially independently propagating deconfined 
spinons, each carrying spin 1/2. The 2-spinon contin¬ 
uum is maximally broad at = tt and shrinks to zero 
at k = 0. This leads to larger contributions to the time 
correlations from the continuum close to fc = tt, which is 
directly visible in the QMC data for long chains as, e.g.. 
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a slower rate of convergence of Rk{t) to a constant. For 
example, i? 27 r/z ,(0 converges much faster to its asymp¬ 
totic constant value compared to i? 7 r(t) for a chain of 
size L = 64, as shown in Fig. [21 In the inset of the same 
figure we also show results for k = tt + 2tt/L, which is 
the momentum we use to extract c as discussed further 
below. Here the chain is longer, L = 200, and Rk{t) does 
not converge sufhciently before the error bars become too 
large. However, the method of fitting a simplified spec¬ 
tral function to QT^it) still works very well and delivers a 
gap consistent with Rk(t), but with much better statis¬ 
tical precision. We therefore exclusively uses the fitting 
method to obtain the results to be discussed next. 

In a finite chain, the ground state is non-degenerate 
and has momentum k = 0. However, there is also a quasi¬ 
degenerate state with k = tt, which is obtained by adding 
an Umklapp to the true ground state, and this state be¬ 
comes exactly degenerate with the k = 0 ground state 
when L — ooi^ There is also a corresponding finite-size 
shift in the excitations close to fc = tt, which is character¬ 
istic of the Heisenberg chain but not present in the model 
in higher dimensions. The lowest triplet gap in the neigh¬ 
borhood of k = TT, where we will use only k = tt + 2tt/L, 
behaves for large L as^ 

O-TT 

A(7r -I- 27r/L) = A(7r) -I- c(T) —, (18) 

Lj 

where A(7r) ~ 1/L but with a multiplicative logarithmic 
correction. In Fig |3] we show the behavior of the corre¬ 
sponding velocity estimate, 

c(L) = ^[A(7r-k27r/L)-A(7r)], (19) 

as a function of the inverse chain length. A smooth mono¬ 
tonic (asymptotically linear) approach to the known ve¬ 
locity, c = 7r/2 can be observed as L —^ oo. There are 
no signs of any remaining log corrections as a function of 
1/L in this estimate. Our results are consistent with a 
remaining linear finite-size correction. 

The velocity of the spinons can also be extracted from 
the lowest triplet gap at A: = 27r/L by using the simple 
estimator LA(27r/L)/(27r) —?► c as L —>■ oo. Results are 
shown in Fig 0) Also in this case we observe the estimate 
approaching the correct value of c as L —>■ oo, but with 
a non-monotonic behavior with a maximum for L Ri 40 
before an apparently linear asymptotic approach to the 
correct value. It should be noted that the spectral weight 
(before normalizing the correlation function) vanishes as 
A: —>■ 0, which implies that the k = 27r/L correlation 
function computed in the QMC simulations becomes very 
noisy for large system sizes. It is worth noting here that, 
because of the non-monotonic behavior seen in Fig 01 a 
calculation using exact (numerical) diagonalization of the 
Hamiltonian would appear inconsistent with the known 
velocity, because large enough system sizes required to 
go beyond the maximum at L 40 cannot be reached. 



FIG. 3: (Color online) Finite-size velocity estimates for the 
Heisenberg spin chain obtained from triplet gaps at and in 
the neighborhood of k = tt according to Eq. m, shown as a 
function of the inverse system size. The velocity approaches 
the known value c = 7r/2 (indicated by the horizontal line) 
for large systems. The line through the larger system sizes is 
a guide to the eye. 



FIG. 4: (Color online) Velocity estimate for the Heisenberg 
chain from the triplet gap close to A: = 0; c{L) = A{k)/k, 
where k — 2tt/L, shown as a function of the inverse system 
size. The known velocity c = 7r/2 in the thermodynamic 
limit is indicated by the horizontal line. The line through 
the larger system sizes is not a fit but is drawn to match the 
known velocity and the data for those system sizes. 

B. Exponential-form projection 

In this section we will discuss a different way of ex¬ 
tracting the gap, using a systematic way to analyze mo¬ 
ments of the spectrum based on information contained in 
the imaginary-time correlations. Such a scheme was re¬ 
cently introduced within T > 0 QMC calculations;^ and 
we here present a generalization to projector QMC. In 
addition, we discuss improvements in the extrapolations 
required to obtain unbiased results. 

To make the connection with the previous version 
of the method more transparent we will here use the 
exponential-form projection with continuous imaginary 
time. The ground state is again projected out of a trial 
state in the valence-bond basis, but now using 
instead of (—H)™. Taylor expanding the exponential 
(as was also done previously in projector QMC, e.g., in 
Ref. [ 40 I 1 . the scheme then closely resembles T > 0 SSE 






QMC algorithmsi^ The normalization, which replaces 
the partition function in T > 0 methods, is expressed as 


Z' = 


= {tjjt\Trexp J dri/(r)^ \tjjt) 

= (^*|1+ ^(-ir / dn / dT 2 


( 20 ) 


/ dTmY\_H{Ti)\'ll}t) 


where indicates time ordering and H{t) is a Hamil¬ 
tonian acting at imaginary time r. If the Hamiltonian 
H{t) here is actually time dependent, we obtain the non- 
equilibrium QMC scheme developed in Ref. |4l|. Here 
H is time-independent, and we formally use the time- 
dependent formalism only as a convenient way of access¬ 
ing imaginary time in the configuration space directly. 
In principle, we can also use the time dependence cor¬ 
responding to the interaction representation^^— where 
only the x- and y— parts of the interaction appears in 
the operator product, but here we stay close to the for¬ 
mulation also used in the SSE method and in the H- 
power method discussed above, and expand with the full 
Hamiltonian. Note that the iJ-power method is exactly 
recovered if the time integrals are completed, and one can 
in practice also equivalently use the power method and 
just generate the ordered time sequences at random;^ 
or generate them in the process of the updates as in the 
world-line continuous time algorithm— 


1. Generalized Moment Method 


Our approach introduced here to extract the gap is a 
generalization to projection QMC of the moment method 
proposed recently— The Fourier-transformed dynamical 
correlation is exploited to derive a sequence of gap es¬ 
timators. We will first explain how the moment of the 
dynamical correlation function is related to the gap, and 
then derive increasingly precise gap estimators. Finally, 
we will show how these estimates are extrapolated to the 
limit in which the exact gap is recovered. 

We will use the following dynamical correlation func¬ 
tion computed with the exponential projector: 


C(t,tq) 




e>i 


where we use the definitions 


\p,q 

= CpCq{p\A^\i){e\A\q), 

(22) 

be 

= befip = co(^ A 0)p, 

(23) 

A^ 

II 

1 

(24) 


We assume that H|0) jP 0, (0|H|0) = 0, and Ai > 0, all 
of which apply here. 

Let us first consider the moment of the asymptotic 
correlation function in the limit f3,TQ —> oo: 


JOC _ 


noo 

/ dTT^Cir^ro) 

Jo 

nOO 


= Y. 




^>1 

n\ 


^n+l 


nl {n » 1). 


(25) 


Then the lowest gap can, in principle, be obtained using 
an appropriate ratio, e.g., 




i: 


Ai 


oo). 


n+1 


(26) 


In practice, however, the projection time (3 cannot be in¬ 
finite in simulations, and we have to consider carefully 
the effects of finite /3. In addition, the range of integra¬ 
tion, T G (0,/3int), is different from (less than) /3 because 
the imaginary-time correlations are measured from the 
reference point tq = (3/2 at the center of the projected 
trial state [noting that Fq. (1^ corresponds to projec¬ 
tion of the bra and ket state with e“^/^], or, alterna¬ 
tively, between two points located symmetrically within 
the projection range r G (0,/3int), and one has to stay 
well away from the boundaries (trial states) for unbiased 
measurements. Moreover, in many cases the statistical 
errors grow too large at large times, as mentioned previ¬ 
ously. Thus, in practice /3int < (3/2. 

The moments for finite (3 and /3int < /3/2 take the form 


L. = 


drr”C(T, To) 




1 rPint 

- / 


£,p 


m 


m—O 


where 


di,p = e 


A^^p — Ei^ Epj 


(27) 


(28) 

(29) 


{13, To -t oo), 


( 21 ) 
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and 


P{m) 


ml 


is a properly normalized Poisson distribution; 


m—0 


(30) 


(31) 


Owing to the finite integration range, the Poisson term 
does not vanish. As a consequence, the ratio of the mo¬ 
ments for finite /3 does not contain information of the 
gap in the limit n —>■ oo. Instead, we have a completely 
different limiting behavior, 


72 “t" 2 1 

^n+1 /5int(^“t“ 1) Pint 


(n —>■ oo). 


(32) 


independent of the gap. We can overcome this difficulty 
and devise a proper gap estimator by using the Fourier 
transformation. Here we express the moments in a dif¬ 
ferent way using the following expansion; 


^2n 

(2n)! 


TOO 

^2n-l 

( 2 n- 1 )! 


lim 

ftnt.T0->OO 


2n ^ ^ ^n,fc,0-^(^fc) 
fc =0 


lim 

/3int.T0->OO 




n 

fc=l 


J{u}k) 

UJk 


(33) 

(34) 


where R{ujk) and J{ujk) are the real and imaginary parts, 
respectively, of the Fourier-transformed correlation func¬ 
tion, i.e.. 


dr 6*™*= (^(t. To) = R{ujk) +*'/(wfc). 


(35) 


where uik = 27rfc//3int {k G Z), and the key coefficients 
are 0 : 14,1 = 1 and 



FIG. 5: (Color online) Extrapolation of A(n —>■ 00 , dint) for 
dint = 6 (diamonds), 7.68 (lower triangles), 9.6 (upper trian¬ 
gles), 19.2 (squares), and 32 (circles). A quadratic polynomial 
in 1/n was fitted to the data points for each dint and the ex¬ 
trapolated n ^ 00 value with error bar is shown in each case. 


Remarkably, this estimator is asymptotically unbiased 
and the limits are interchangeable (see Appendix for the 
analytical derivation): 


Ai 


lim lim 

n-)-oo /3i„t,To->oo 


A 




lim lim A 

Sinti’7‘0—>-oo n—>oo 


(?T.,^int) * 


(38) 


Due to the commuting limits, observing convergence of 
the estimator (1571) in large n and dint taken in any conve¬ 
nient fashion will deliver an unbiased result for the gap. 

In principle there is also a dependence on the d value 
used in the exponential projection, in addition to the 
dependence on n and dint- Above we have assumed that d 
is sufficiently large for quantities computed by integration 
up to dint have converged to the d 00 limit, and in 
the QMC simulations we also monitor this convergence. 
Naturally, the value of d required grows with dint- 


2. Performance Check for the Heisenberg Chain 


^n,k,m — r - (56) 

n (^+.7')(fc--7') 

m<j<n,j^k 

When deriving these equations we have considered the 
expansion in ojk on the right-hand side of Eq. (1331) and 
Eq. (IMl) . where the lowest orders then cancel. The coef¬ 
ficients Xn,k,m can be solved for by using the inverse of 
the Vandermonde matrix.— 

Combining the results above, we obtain the following 
improved gap estimator: 


We have checked the validity of the estimator (1571) for 
the Heisenberg chain of 16 spins with periodic boundary 
conditions. The projection length d is set large enough 
to see the asymptotic behavior of the dynamic correla¬ 
tions on the imaginary-time axis. The correlations were 
measured from the center of the QMC configuration; that 
is, from tq = d/2- As in the previous section, the loop 
algorithm and improved correlation-function estimator- 
are used. We consider the triplet gap at A: = tt and 
the operator A in Eq. m is then explicitly given by 
A = J2r Sr(-I)’’ and 


E 




A 2 k—1 

^(n,/3int) — n 


J (^/c) 

^k 






^ (39) 

OL—x,y,z r,r' 

The Fourier-transformed correlation functions (1351) were 
directly measured in the simulation, and no discretiza¬ 
tion error is introduced. The gap estimators (|57)) for 
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FIG. 6: (Color online) Extrapolation of the gap for Pint ^ oo. 
An exponential function /(/3int) = A + a x exp(— 6/3int) (curve) 
has been fitted to the data points, where a and b are positive 
real numbers. The exact gap value calculated from the full 
diagonalization is shown by the horizontal line. 

several n and Pint were calculated and the errors esti¬ 
mated by the jackknife analysis^ (bootstrapping would 
produce statistically equivalent result). 

As shown above, our estimator (l37l) converges to the 
exact gap in the limit of infinite n and /3int- In prac¬ 
tice, if good convergence is observed within statistical 
errors, a gap value obtained from a sufficiently large n 
and Pint within the converged range can be used as the 
final estimation, or some appropriate functional form can 
be used for extrapolation. However, an important issue 
here is that the statistical errors grow with the two pa¬ 
rameters. We therefore inevitably encounter a trade-off 
problem between systematical and statistical extrapola¬ 
tion errors. Here we will demonstrate that the statistical 
error can be optimized by proper extrapolations while 
keeping the estimation unbiased. In the procedure used 
below, the extrapolation is taken for n —> oo first and 
then for Pint —>■ oo, though, as discussed above, other 
ways to accomplish the limit are also possible. 

Examples, based on (~ 10®) Monte Carlo samples, 
of n —> oo extrapolations are shown in Fig. [5j The lead¬ 
ing finite-n correction is linear in 1/n [in accord with 
Eq. (IA4I) and (lAlOl) in the Appendix], and we include 
also a quadratic term to obtain good fits to the data for 
the full n-range. The limit Pint —> oo is taken next us¬ 
ing the resulting values of the n —> oo extrapolations, 
as shown in Fig [51 Here an exponential function is used 
for the data fit (again, according to results derived in 
the Appendix) to extrapolate the final result for the gap. 
Though other ways of extrapolating to n. Pint oo are 
possible, the above protocol is convenient because the ex¬ 
trapolation for n, which is taken first, is relatively easier 
than that for Pint- 

As a test of the unbiased nature of the extrapolation 
scheme, distributions of the relative gap error are shown 
in Fig. [71 Histograms were collected based on 2048 in¬ 
dependent simulations of the L = 16 Heisenberg chain. 
Results based on the extrapolation procedure discussed 



Relative error 

FIG. 7: Histograms of the density of the relative error of 
the (n —^ oo. Pint —>■ oo) extrapolated gap, along with those 
based on the estimator (I37II for fixed (n, Pint) = (1)32) and 
(10,32). The error was calculated as (A — A)/A, where A 
is the generalized-moment estimated value from 2^® (~ 10®) 
Monte Carlo samples and A is the exact value for the L = 
16 Heisenberg chain. The histograms include results of 2048 
independent simulations. 

above are compared with those of individual gap estima¬ 
tors for (n,/3int) = (1,32) and (10,32). The (1,32) esti¬ 
mator clearly has a non-zero systematic error remaining, 
which is similar to (but smaller than) the conventional 
second moment estimatorThe (10,32) estimator has 
a small enough systematical error (the histogram being 
centered very close to zero) but has a large statistical er¬ 
ror (wide distribution). The extrapolated estimation is 
unbiased and has a small statistical error. 


III. THE VELOCITY FROM WINDING 
NUMBER FLUCTUATIONS 

We here discuss how to use winding numbers to com¬ 
pute the velocity. This method has been known for some 
tiinefi^ and was recently applied to the 2D XY modeli^ 
and the 2D Heisenberg model,^ We begin by briefly recol¬ 
lecting some key aspects about winding numbers in finite- 
temperature QMC simulations, in particular the SSE 
method we use for these calculations. We then present 
results for the Heisenberg chain. 

A. Winding numbers in QMC simulations 

QMC simulations at finite temperature are based on 
some mapping of the partition function of a quantum sys¬ 
tem in D dimensions to an effectively equivalent {D + 1)- 
dimensional classical system, where the additional di¬ 
mension of the configurations corresponds to imaginary 
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(Euclidean) time. The effective length of the system in 
the time dimension is c/3, where /3 = 1/T (setting ks = 1) 
and the configurations are time-periodic. For a system 
with conserved particle number, which in the case of spin 
model corresponds to conserved magnetization along the 
quantization (z) axis, imposing periodic boundary in the 
spatial dimensions leads to another, topological number 
associated with the configurations—the winding number 
representing permutations when particles are propagated 
once or multiple times around the periodic system. The 
winding numbers were first used in QMC calculations of 
the superfiuid stiffness of interacting bosons 

The SSE QMC algorithm j^i^^d^ is based on a Taylor 
expansion of the imaginary-time evolution operator (the 
Boltzmann operator) , and similarities with the pro¬ 
jector approach discussed in the previous section were al¬ 
ready pointed out. Each SSE configuration is associated 
with some power i?" of the Hamiltonian propagating a 
basis state (here in the standard computational basis of z 
spin components), and these powers are sampled stochas¬ 
tically to all contributing orders. The trace over all basis 
state is also sampled. The average expansion power (n) 
in this procedure is proportional to (3; (n) = /3{H) oc PN. 
In simulations, the state propagation is broken up into 
individual paths corresponding to strings of n of the in¬ 
dividual local terms of the Hamiltonian, forming succes¬ 
sions of n evolving basis states, similar to those in path 
integrals. For a Heisenberg model the Hamiltonian terms 
are the diagonal operators SfSj (in practice with a con¬ 
stant added) and off-diagonal S~ + S~ operators, 
the latter of which transport spin and are associated with 
currents Ja = ±1 in the lattice direction a corresponding 
to the site-pair i,j. The winding number in the a lattice 
direction is defined in terms of the currents as 

1 " 

W^a = ^5^Ja(p), (40) 

where the index p corresponds to the location of the 
transport “event” in the string of n operators and La 
is the length of the lattice in the a direction. Defined in 
this way, the winding numbers are integers. It should be 
noted that the definition of the winding number is ex¬ 
actly the same in SSE and world-line methods^i^ and 
one can also think of the SSE configurations as consist¬ 
ing of world lines (for up and down spins in the case of 
S = 1/2 quantum spin systems). 

The spatial winding number Wr measures the net spin 
transported around the periodic lattice in the r direction 
in the course of the periodic propagation in imaginary 
time. Equivalently, this is the number of world lines (up 
ones minus down ones divided by two) crossing through 
a plane drawn along the time axis perpendicular to the 
r axis. Since the total z magnetization is conserved, 
one can also think of as a winding number; the net 
number of world lines crossing a plane drawn at an ar¬ 
bitrary time point perpendicular to the time axis, which 
is just the magnetization computed in the stored basis 


state; 

N 

Wr = M, = J2s/. (41) 

i=l 

The expectation values of the squared winding numbers 
(i.e., the winding number fluctuations) are related to two 
important thermodynamic quantities; the spin stiffness, 

Ps = ^{{W^) + {WS)), ( 42 ) 

and the uniform magnetic susceptibility, 

X = V <"?> = V 

The technicalities of implementing these observables in 
SSE calculations have been discussed extensively in the 
literature; see Ref. 113 for a recent review. 

B. The cubic criterion and the velocity 

In the high-temperature limit T —^ oo, the magnetiza¬ 
tion fluctuations of any system are maximized and there¬ 
fore (W/) > 0 according to Eq. (l43l) . For an unfrustrated 
antiferromagnet the ground state is a singlet, and, on ac¬ 
count of the presence of a singlet-triplet finite-size gap, 
{W/} —>■ 0 when T —>■ 0 for any finite system. In contrast 
to these limits of (W^), for the spatial winding number 
in direction r (r = x,y,...) we have (W//) —>■ 0 when 
T —>■ oo on account of there being no quantum fluctua¬ 
tions when the imaginary-time length /3 —>■ 0 and there 
are no contributions from expansion powers n > 0 (in 
the case of SSE—in world-line methods there will simi¬ 
larly be no transport events causing shifts of the world 
lines). In the limit T —0, for a system with long-range 
order (or a ’’quasi-ordered” ID system with power-law 
decaying correlations), the stiffness constant converges 
to a non-zero value for any L, and according to Eq. 
we must then have a divergence {W/) ^ 1/T. These 
different behaviors of the spatial and temporal winding 
numbers versus temperature guarantees that there is a 
crossing point, 

{w/{n) = {w^{n), (44) 

at some unique value of /3 = 13* {L) for given size L. 

The winding numbers characterize global fluctuations 
of the system in the different spatial and temporal direc¬ 
tions. It is then natural to define a system as having 
cubic space-time geometry when Eq. (ITil) holds (with 
the lattice length L the same in all spatial directions). 
The aspect ratio L/(3* should then directly correspond to 
the velocity of the long-wave-length excitations. In some 
cases this can be shown directly based on low-energy field 
theor y^d^d^ even in the absence of such descriptions 
the arguments are very general and one can expect the 
conclusion c = L/f3* to always hold for a system with 
linear dispersion, though we are not aware of any formal 
proofs in the general case. 
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FIG. 8: (Color online) Difference between the spatial and 
temporal winding numbers versus the inverse temperature in 
simnlations of a 64-site Heisenberg chain. A second-order 
polynomial has been fitted to the data points, and this curve 
is used to determine the value /3 = /3* at which cubic condi¬ 
tion (Mf?) = (yVr) holds (i.e., intersection with the horizontal 
dashed line). 



FIG. 9: (Color online) Size dependence of the velocity ob¬ 
tained with the cubic-criterion for the Heisenberg chain, along 
with a fit corresponding to a leading ~ l/L^ size correction. 
The dashed line is at the rigorously known velocity c = 7r/2. 


C. Test on the Heisenberg chain 


To find the point where the cubic criterion (W^) = 
{W"^) is satisfied, we simulate a system at several values 
of (3 in the region where (W^) ~ (Wt) based on initial 
explorations and knowledge of the approximate value of 
the velocity. We fit a low-order polynomial (typically 
second- or cubic-order) to the difference (W)^) —(FF^) and 
solve the resulting equation for the /3-value for which the 
cubic criterion is satisfied. This procedure is illustrated 
in Fig. |8] for a Heisenberg chain with L = 64 spins. From 
this procedure we obtain c{L) = L/j3*{L), which can be 
extrapolated to L —?► oo. With independent data points, 
the statistical error of c{L) at fixed L can be determined 
by repeating the procedure multiple times with added 
Gaussian noise of standard deviation equal to the error 
bars of the data points, whence the standard deviation 
of the extracted crossing point is the error bar. 


Fig. [9] shows results of such a procedure for several 
chain lengths L, graphed versus 1/L. We are not aware 
of any theoretical predictions for the size dependence of 
this definition of c(L), but the data for the larger systems 
are well described by a constant (the final infinite-size 
value of c) plus a term proportional to 1/L^. Using this 
fitting form leads to a value of c completely consistent 
with the known value c = 7r/2, as is clear from Fig. [HI 


IV. HYDRODYNAMIC RELATIONSHIP 

A well known way to extract the velocity in an anti- 
ferromagnet is to use the analogue of a hydrodynamic 
relationship between the velocity, the spin stiffness (he- 
licity modulus), and the transversal susceptibility (which 
is the analogue of a mass density) 



In a QMC calculation in which the spin-rotation sym¬ 
metry is not explicitly broken, one cannot compute x_l 
directly, but one can use the fact that the rotationally- 
averaged uniform susceptibility x computed in the z ba¬ 
sis, as in Eq. (PI). becomes 2/3 of a transversal (e.g., 
x) component when T —> 0 in the thermodynamic limit 
(since the longitudinal component vanishes at T = 0). 
Thus, one can obtain xx (3/2)x by taking the limit 
L —>■ oo before the T —0 limit, while taking the limits in 
the opposite order does not work because then x —0 due 
to the finite-size gap between the = 0 and > 0 
magnetization sectors. In contrast, for ps the limit T —>■ 0 
has to be taken before L —>■ oo, because the system in the 
thermodynamic limit only has stiffness (Neel order) ex¬ 
actly at T = 0. In this case, too, a factor 3/2 has to be 
included in the finite-size estimate to account for rota¬ 
tional averaging of the relevant transversal components. 

The procedure to obtain ps in the thermodynamic limit 
is relatively straight-forward with an extrapolation of 
Ps{L,T —>■ 0) using a polynomial in 1/L, while the ex¬ 
trapolations requiring L —>■ oo first in the x calculation is 
more cumbersome. Results obtained for c in this way^ 
have large error bars compared to the result of the wind¬ 
ing number method presented above in Sec. IV Al 

In order to define a finite-size velocity estimate c{L) 
based on Eq. (I45|) directly in the T = 0 limit we here use 
a modification of the relationship. We use the suscepti¬ 
bility at finite momentum q, 

Xin) = ^ (46) 

where the magnetization at non-zero momentum is the 
same as A in Eq. 0 but without the normalization by 
We can use the smallest momentum q = 211 jL 
for a given system size and approach the q = 0 limit as 
L —> oo. If we here define Ps{L) without the rotational 
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FIG. 10: (Color online) Size dependence of the velocity of 
excitations in the Heisenberg chain, defined according to the 
hydrodynamic relationship in the form m- The apparent 
poor convergence to the exact result c = 7r/2 when L —>■ oo 
can be explained by logarithmic corrections. The maximum 
can be related to the anomaly in the excitation spectrum at 
k = 2ir/L, which causes the maximum at the same L in Fig. [4] 

factor 3/2 discussed above, then the same rotational fac¬ 
tor should also not be used in the susceptibility, and we 
define the velocity for finite size as 

= \^x(/=27!/L)- 

We will compute this quantity using SSE simulations at 
sufficiently large j3 to achieve ground state convergence 
for each L studied. 


A. Test on the Heisenberg chain 

The hydrodynamic relationship (1451) is normally ap¬ 
plied in the magnetically ordered Neel state and one may 
then question its use in a critical system. In the case 
of the Heisenberg chain we can rigorously see that it is 
a valid way to extract c. Since there is no long-range 
order, there is no distinction between longitudinal and 
transversal modes, but Eq. (1451) defined with rotation- 
ally averaged quantities should remain valid. According 
to the exact Bethe Ansatz solution, in the thermody¬ 
namic limit we have, in the units with J = 1 and all 
relevant constants set to 1, ps = 1/4 (derived in Ref. nil) 
and x ((7 —>■ 0) = I/tt^ (see., e.g.. Ref. Is^l and, thus, ac¬ 
cording to Eq. dUD, we obtain the correct result c = 7r/2. 

Eor hnite system size, it is well known that both x 
and Pa exhibit strong logarithmic corrections to their 
infinite-size values which can be traced to the presence of 
marginal operators and it is very difficult to extrapolate 
them; see, e.g.. Ref. [H for ps and Ref. [H for x (where in 
the latter case the logarithmic correction in temperature 
is discussed but one should expect similar corrections for 
g —>■ 0 at T = 0). The logarithmic corrections do not 
cancel in the ratio of the two quantities in Eq. (SB, and 


as a consequence we also find that it is difficult to extrap¬ 
olate c{L) precisely to the thermodynamic limit. Results 
are shown in Fig. (TU] Although it is not possible to fully 
extrapolate to the exact results based on these data for 
system size up to A = 512, the result for the largest size 
nevertheless deviates by only 0 . 2 % from the exact result. 

It is interesting to note that also this definition of c(A) 
exhibits a non-monotonic behavior, with a maximum at 
L Ki 40, the same as in the c{L) value obtained previ¬ 
ously from the gap at fc = 27r/L in Fig. 21 In the latter 
case, the maximum corresponds to an elevated excitation 
energy at fc = 27r/A, which one should expect to lead to 
a reduction in x(27r/i) because this quantity is given by 
a sum rule of S{k,u})/ui. By this sum-rule, in a single¬ 
mode approximation a larger uik would lead to a smaller 
x(fc), and even beyond the single-mode approximation 
one should expect such an effect because the dominant 
spectral weight is at the gap. If pa is not appreciably 
affected by this finite-size anomaly, then indeed an ele¬ 
vated gap and reduced x(27r/A) in Eq. (ITT)) can explain 
the maximum in Fig. HOI 

In light of the logarithmic corrections to c{L) defined 
according to Eq. (251) , the apparent complete lack of any 
challenging corrections or non-monotonic behavior in the 
definition based on the cubic criterion for winding num¬ 
bers (Fig. [S]) is even more remarkable, since also this 
estimate involves quantities directly related to the sus¬ 
ceptibility (temporal winding number) and spin stiffness 
(spatial winding number). There are also no apparent 
logarithmic corrections in the velocities defined based on 
gaps in Figs. [3] and 21 In Fig. 21 which is based on the 
triplet gap close to tt, the logarithmic corrections are 
avoided by subtracting the k = tt gap in Eq. o, while 
the gap close to fc = 0 used in Fig. 21 is not expected to 
be affected by logarithms. The latter gap at fc = jL 
scales as ck + b/V^ for large L according to our results in 
Fig. 21 which can be explained from the known disper¬ 
sion Afc = c sin(fc) = ck-\- dk^ -I-..., if there is a finite-size 
correction ~ 1/L^ (due to irrelevant fields only). 


V. 2D HEISENBERG MODEL 

The 2D spin-1/2 Heisenberg model on the square lat¬ 
tice spontaneously breaks spin rotation symmetry at 
T = 0 in the thermodynamic limit. This leads to gapless 
linearly dispersing Goldstone modes in the vicinity of the 
momenta (0,0) and (tt,tt). The ground state in a finite 
periodic system is a total spin singlet. The long-range 
antiferromagnetic order in the thermodynamic limit is 
reflected in the energies of the S' > 0 quantum ro¬ 
tor states^ which collapse onto the ground state as 
A 5 ~ I/LP' much faster than than the spin wave exci¬ 
tations, A ^ 1/A. The rotor states, thus, become degen¬ 
erate with the ground state as the system size increases, 
and have momenta k = ( 0 , 0 ) and (tt, tt) for even and odd 
S respectively. Combinations of rotor states with S up 
to ~ A can then be formed which are ground states with 
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FIG. 11: (Color online) Size dependence of the velocity de¬ 
fined using the cubic-criterion for the 2D Heisenberg model, 
along with an fifth-order polynomial fit (including all system 
sizes L > 6) without linear and quadratic terms. The extrap¬ 
olated velocity is c = 1.65847(4). The data for the largest 
system sizes are shown on a more detailed scale in the inset. 
The goodness of the fit is x^/dof « 0.8. 

fixed direction of the Neel vector (the staggered mag¬ 
netization), thus allowing for symmetry breaking in the 
thermodynamic limit. Here we provide an accurate de¬ 
termination of the spin-wave velocity and also directly 
investigate the scaling of the rotor state. 


A. Velocity from winding nnmbers 

Since we expect the winding numbers with the cubic 
criterion to provide the best determination of c we start 
with this approach. Results based on the procedures dis¬ 
cussed in the preceding section are shown in Fig. [TT] ver¬ 
sus the inverse system length. Carrying out polynomial 
fits, we find that no linear and quadratic terms are re¬ 
quired. Fourth- and fifth-order polynomials fits excluding 
these terms and using all the L > 6 data give almost iden¬ 
tical results for the L —^ oo extrapolated velocity, with 
only the error bar somewhat larger for the fifth-order fit. 
The hgure shows the hfth-order fit and the extrapolated 
velocity with it is c = 1.65847(4). While we do not know 
the physical reason for the leading cubic correction, we 
use it as the simplest empirical description of the data. 
Including also a quadratic term, it comes out equal to 0 
within statistical errors and the extrapolated result does 
not change appreciably. 

The above value of c agrees within errors bars with the 
recent result using the same method by Jiang and Wiese^ 
but our error bar is almost an order of magnitude smaller. 
We note that in Ref. @ no systematic extrapolation was 
carried out to the thermodynamic limit—instead an av¬ 
erage was taken of results for system sizes in the range 
L G [24,64]. Looking at the data in Fig. [TT]it is clear 
that, with the small error bars on the SSE data we have 
achieved here, an extrapolation is necessary to obtain 
a result with no remaining finite-size effects. To our 



FIG. 12: (Color online) Size dependence of the lowest triplet 
gap A(7r, tt) multiplied by for the 2D Heisenberg model. 
This analysis shows that the gap scales as 1/L^ for large L 
and an extrapolation based on a polynomial in 1/L gives the 
estimate l/2x± = 7.62(2) for the uniform susceptibility. 

knowledge, the above result is the most precise spin-wave 
velocity reported to date for the 2D Heisenberg model. 
Spin-wave theory with corrections up to order 1/S'^ gave 
c = 1.6638(3))^ where the uncertainty 3 in the last digit 
reflects estimated numerical errors from evaluations of 
challenging integrals. Thus, to this order, the spin-wave 
result deviates by only 0.3% from the correct value. 

B. Gap scaling of rotor states 

The quantum rotor excitation gap can be directly ac¬ 
cessed by measuring the lowest triplet gap at k = (tt, tt). 
This energy scale is related to the uniform (transverse) 
magnetic susceptibility x_l as)^ 

= ( 48 ) 

From the behavior of the lowest triplet gap up to sys¬ 
tem sizes L = 64, as shown in Fig. [1^ we indeed ob¬ 
serve that A(7r,7r) ~ 1/L^ at large L, but there are 
also large corrections which we fit with additional higher- 
order powers of 1/L. The extrapolation to infinite size 
gives I/2 xj_ = 7.62(2). This is consistent with the 
value of susceptibility obtained using QMC calculations 
in small external magnetic fields to extract gaps between 
different spin sectors 

C. Velocity from gaps 

The velocity of the spin waves can be estimated by 
measuring the triplet gap in the vicinity of (7r,7r) and 
(0,0). We here choose to measure the triplet gap at 
ki = (tt + 27r/L, tt) [or, equivalently, at (tt, tt + 27r/L), 
which we use for averaging] since the lowest triplet ex¬ 
citation energy is at (tt, tt) and ki is the closest allowed 
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FIG. 13: (Color online) Extraction of the gap at fc = 
{it -\-2'k!L,!^) based on power-projected correlation functions, 
shown for the largest 2D Heisenberg system studied; L = 64. 

(a) The ratio Rk{t), Eq. (IIOII converges to a hnite gap value 
but the error bars are large for long imaginary times. Fitting 
Qk{t) to four S functions according to Eq. (Hill , as shown in 

(b) , provides a more reliable gap estimation, as illustrated in 
the inset of (a). The thickness of the solid line here approx¬ 
imately represents the error bar. The inset of (b) shows the 
distribution of gap values obtained from a large number of 
different bootstrap samples of the QMC data, from which the 
error bar of the gap estimate was computed. 

wave-vector to (tt, tt) for a periodic system with linear 
size L. For this model we have carried out only i7-power 
QMC simulations. We illustrate the gap extraction with 
both the ratio method and the simplified spectral func¬ 
tion with fitted (5-functions in Fig. [T3l using the largest 
system sizes considered; L = 64. Again, working with 
the spectral function produces much more stable results, 
though clearly the correlation ratio Rk also converges to 
a constant consistent with the same gap. 

Since in the thermodynamic limit, the excitation en¬ 
ergy of the spin waves equals E(k) = c|k — (tt, 7r)| in the 
vicinity of (7r,7r), where c is the spin-wave velocity, the 
estimator 

c(L) = — A(7r -b 27r/L, tt) (49) 

ZTT 

should converge to c as L —> oo. We graph this quantity 
in Fig. [TU and it converges to the value of c obtained 
above from the winding-number method when L —>■ cx). 



FIG. 14: (Color online) The spin-wave velocity estimator c{L) 
defined in Eq. (1491) shown for the 2D Heisenberg model as a 
function of 1/L. The velocity for large L is consistent with the 
value obtained from the winding numbers ('Fig. 1111) . shown as 
a horizontal line. The line through the values for the larger 
system sizes is a guide to the eye. 

Note that there is again a non-monotonicity as a function 
of L, similar to the case of the Heisenberg chain in Fig.|4l 
In the latter case this behavior only was observed in the 
gap extracted close to = 0, however, while here we have 
used the spectrum close to (7r,7r). 

D. Hydrodynamic relationship 

We next consider the definition of c{L) by the hydrody¬ 
namic relationship, Eq. gZl). We again go to sufficiently 
low temperature in SSE simulations for any remaining 
finite-T corrections to be negligible compared to the sta¬ 
tistical errors, which typically meant /3 = 8L. Results 
are graphed in Fig. 1151 Here again we expect the finite 
size corrections in the Neel state to be described asymp¬ 
totically by a polynomial in 1/L, but a rather high or¬ 
der of the polynomial is required to fit the data for the 
smaller system sizes. The behavior for the largest sizes 
again indicate (as in the case of the winding-number cal¬ 
culation discussed above) that the leading correction in 
1/L is cubic. To get a statistically acceptable fit for all 
L > 8 data we include terms up to order 6, which gives 
c = 1.65875(10), which deviates by approximately 2.5 er¬ 
ror bars from the statistically much more precise value 
obtained in Sec. IV Al based on the winding-number fluctu¬ 
ations. By comparing Figs. [TT] and fTKl it is clear that the 
analysis in the former is more reliable, with a larger num¬ 
ber of data points used and an overall much weaker size 
dependence. The above result from the hydrodynamic 
relationship is still likely somewhat affected by system¬ 
atical errors, as the behavior still appears to be flatten¬ 
ing out more than the fit suggests, and the fit even in¬ 
cluding the 6th-order term is only marginally good, with 
X^/dof ~ 1.7. It would clearly be desirable to go to larger 
system sizes, but given the much better behavior of the 
winding number data it is already clear that this is the 
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FIG. 15: (Color online) Finite-size spin-wave velocity defined 
using the modified hydrodynamic relationship (I47II . The ex¬ 
trapolation to infinite size is done using a sixth-order polyno¬ 
mial in l/L without linear and quadratic terms (solid curve). 
The inset shows the data for the largest systems on a more 
detailed scale. The point (blue) with error bar close to the 
y-axis indicates the statistical error of the L —>■ oo extrapo¬ 
lation. The barely separated horizontal lines shows the value 
of c plus-minus one error bar obtained in Sec. IV Al using the 
winding number method. 

preferred method. 

VI. BILAYER HEISENBERG MODEL 

We now consider the S' = 1/2 bilayer Heisenberg model 
which is a prototypical system to study quantum phase 
transitions in 2D.— The Hamiltonian is given by 

H = JiY. (Sm • S,.i + S ,.2 • S,. 2 ) + J 2 ^ S,.i • S,. 2 , (50) 

{i.j) i=l 

where represents a S = 1/2 spin operator at site i 
of layer a (a = 1,2), and (i,j) denotes a pair of nearest 
neighbor sites on the square lattice oi L x L sites with 
periodic boundary conditions. Both the couplings Ji, J2 
are positive (antiferromagnetic). As the ratio g = J2/J1 
is increased, there is a destruction of long-range Neel or¬ 
der at a critical gc, beyond which the system enters a 
disordered state with no broken symmetries. The best 
value available for the location of the critical point is 
Qc = 2.5220(1)^ and we will use this value below. 

A. Velocity from winding nnmbers 

Results for the velocity based on the cubic winding- 
number criterion at the critical point of the bilayer are 
shown in Fig. 1161 Here we find that a non-integer power- 
law correction describes the data very well for system 
sizes L > 6. A fit of the form c{L) = c+bjL‘^ to the L > 8 
data gives c = 1.9001(2). Note that the fitted curve also 
goes through the L = 6 data point even though this point 



FIG. 16: (Golor online) Size dependence of the velocity com¬ 
puted using the winding-number cubic criterion for the bi¬ 
layer Heisenberg model at the estimated critical coupling ra¬ 
tio g — 2.5220. The curve shows a quadratic fit of the form 
c{L) = c -I- b/L‘^ to the L > 8 data, which gives the extrapo¬ 
lated velocity c = 1.9001(2). 

was not included in the fit. This adds to our confidence of 
this power-law correction. The value of c is a few percent 
smaller than the spin-wave resultiS Csw = 1.96 including 
I/S' corrections at g = 2.51 and Csw ~ 2.03 obtained^ 
to orderl/S^. A more sophisticated treatment beyond 
conventional spin-wave theory, including the effects of 
longitudinal fluctuations close to the critical point, gives 
c = 1.78 when the expression c = 0.705 x g below Eq. (10) 
of Ref. [13 is evaluated with gc = 2.522. While these ana¬ 
lytical values may appear reasonably close to (deviating 
by a few percent from) the presumably correct numeri¬ 
cal value we have obtained here, the deviations are still 
much larger than in the case of the ordered 2D Heisen¬ 
berg model, where the error of spin-wave result is only 
0.3%, as discussed in Sec. IV Al 

The exponent of the correction in the fit in Fig. [TOl is 
a = 1.67(4). It is not clear to us how this exponent relates 
to the standard critical exponents of the 3D 0(3) univer¬ 
sality class of the transition, but the relatively large value 
(larger than 1/v Ki 1.41) suggests that it may involve 
subleading exponents. 

B. Velocity from gaps 

Since the bilayer has two sites per unit cell (a = 1,2 for 
each square-lattice point I), the triplet excitations have 
an extra label in k space which denote the in-phase 
{kz = 0) and out-of-phase {kz = tt) spin excitations of 
the two layers, respectively. In the magnetically ordered 
Neel phase, the triplet excitations are gapless at (0,0,0) 
and (tt, tt, tt), with the lowest excitation being at (tt, tt, tt) 
for a finite system. The spectrum is linear in the neigh¬ 
borhood of both the gapless points, which defines the 
corresponding spin wave velocity c. For a continuous 
phase transition, the spin wave velocity c scales asi^ 

c oc (g - gc) 


iz-l) 


(51) 















17 



FIG. 17: (Color online) The velocity estimator (I52II for the 
Heisenberg bilayer at its critical point shown as a function of 
1/L. Results based on both power-projection and exponen¬ 
tial projection (continuous time) are shown and agree within 
error bars. The velocity estimate obtained from the winding 
numbers (Fig. I16II is shown as the horizontal line. 

where z is the dynamical exponent and v is the correla¬ 
tion length exponent. Thus, for a z = 1 transition as is 
the case for the Heisenberg bilayer a.t g = gc, the velocity 
is regular at the critical point. 

1. Critical point 

For measuring the velocity of the critical modes we 
first study the triplet gap at (tt -|- 27r/L, tt, tt). From the 
linearly of the spectrum at the critical point, we define 
the velocity estimator in analogy with the single-layer 
case as 

c{L) = A{tt + 2tt/ L, n , tt) L /{2tt) . (52) 

The behavior of this quantity as a function of L is shown 
in Fig. [T71 Unlike the case of the Heisenberg chain 
(Fig. [5]) and the single-layer (Fig. ITT)) , the velocity es¬ 
timate for the Heisenberg bilayer at criticality is notably 
higher than winding-numbers estimate, by about 5% as 
L —7^ oo. To extract the gap needed in Eq. (l5^ we 
have used both the power-projection QMC method with 
i5-function fits to the correlation function, as described 
in Sec. Ill HI and the generalized moment method ap¬ 
plied to exponential-projector QMC data, as detailed in 
Sec. IHB II As seen in in Fig. [T71 the two methods give 
results that agree fully within statistical error. 

At first sight, potential reason for the disagreement 
with the winding-number result could be an over¬ 
estimation of the gap extracted from the correlation func¬ 
tion at the critical point. Consider the full dynamic spin 
structure factor defined as 

5k(cc) = ^ |(m|5^(k)|0)|2(5(u; + Eo - Em). (53) 

m 

In the Neel phase the structure factor has a (5-function 
of weight Ad(k) at an energy w(k) which represents the 



FIG. 18: (Color online) Momentum dependent gap, with q 
defined relative to {tt,tt,tt), divided by the respective q. The 
data are shown versus 1/L for q = 27r/4 (open triangles), 
27r/5 (open diamonds), 2tt/6 (open circles), 2tt/8 (solid tri¬ 
angles), 27r/10 (solid diamonds), 27r/12 (solid circles). Linear 
fits in 1/L are shown for each momentum, with only data for 
sufficiently large L included for each q. 

magnon mode, as well as a continuum Ac(k, w) which 
does not extend below a;(k) and which decays rapidly to 
zero as w —>■ oo: 

5'k(w) = Ad(k)(5(w - w(k)) -|- Ac(k, w). (54) 

The velocity derived from a;(k) in the vicinity of (tt, tt, tt) 
is a regular function oi g — g^ and is by definition the cor¬ 
rect velocity at the critical point. However, the weight in 
the ^-function ^^(k) (that represents the magnon mode) 
also smoothly goes to zero as the critical point is ap¬ 
proached, with the spectrum evolving into purely a con¬ 
tinuum reflecting the overdamped critical inagnons. In 
our fitting scheme [see Eqs. m and (ITTl) ] used with the 
power-projection QMC method, where the spectral func¬ 
tion is represented by a small number of S functions, u;(k) 
may be over-estimated, especially for large system sizes, 
when Ad is very small close to the critical point. The 
velocity would then also be over-estimated, as it is in 
Fig. |T7] (assuming that the winding-number result is cor¬ 
rect). A similar distorting effect of the continuum may be 
expected also with the generalized moment method used 
with the exponential-projection QMC data. However, 
the (5-function fits are stable with respect to the number 
of (5-functions used, and the extrapolations used with the 
generalized moment method are also stable. We do not 
see any evidence of remaining effects that could account 
for an overestimation as large as in Fig. [T71 The perfect 
match between the two methods within their statistical 
errors also gives us confidence that the gaps are deter¬ 
mined correctly and the reason for the disagreement with 
the winding-number method must be sought elsewhere. 

To search for a possible flaw in the velocity estimation 
based on Eq. (1521) . we next analyze the details of the 
dispersion relation around (tt, tt, tt), using gaps extracted 
with the generalized moment method. The dispersion 
should be an asymptotically linear function of the wave 
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FIG. 19: (Color online) Extrapolation of the critical bilayer 
velocity using the gaps for q —>■ 0. The infinite-size veloc¬ 
ity is estimated as the limiting value c = 1.899(2), which 
is consistent with the estimation from the winding numbers, 
c = 1.9001(2), indicated by the horizontal line. A linear func¬ 
tion in 1/q^ is used for the extrapolation. The inset shows 
the calculated dispersion. 

number, A(q) = eg, where q is the momentum relative 
to the gapless point (tt, 7r,7r); 

q = k — (tt, TT, tt). (55) 

We should then have limg_>o limi_^oo A(q, L)/g = c. 
There could potentially be an issue with the order of 
the limits g —>■ 0 and L —>■ oo, which with Eq. (l52ll are 
taken simultaneously as we use g = 27r/L. 

To investigate the formally correct limit of taking L —>• 
oo first and then g —>■ 0, we analyze the gaps at fixed 
g and varying L, to find the corresponding gap values 
as a function of the momentum in the thermodynamic 
limit. The L dependence of the gaps at several different 
g-values are shown in Fig. [THl We observe that the size 
correction to the gap for sufficiently large L is linear in 
1 /L and therefore extrapolate the data to infinite size 
using simple line fits. Using these gaps A(g), we finally 
extrapolate the velocity as c = limfc^o A(g)/g, as shown 
in Fig [T21 The behavior is fully consistent with a linear 
dispersion in the limit g —^ 0, and empirically we find that 
the leading correction is of order g^ (i.e., the correction 
to the gap A ~ g is of order g^). The final result for 
the velocity extracted this way is c = 1.899(2), which 
is fully consistent with the result based on the winding 
numbers. Thus, the reason for the previous disagreement 
is indeed that Eq. (l5^ does not represent the correct 
order of limits, and there are no flaws in the extraction 
of the gaps. 

This more detailed analysis also shows the reason for 
the over-estimation of the velocity based on Eq. (l5^ . As 
is clear from Fig [T8j the gaps converge to their infinite 
size values as 1/L for sufficiently large L. We observe 
that the gaps at fixed momentum k, where k is close to 
but not equal to (tt, tt, tt), follows the following behaviour 


for large L: 

A(g, L) = A(g, L -)> oo) -b (gapless) 

A(g,L) = A(g, L —» oo) -b (critical) 

Lj 

A{q,L) = A(g,Loo)-b a(g)e“^^'^^^ (gapped) 

In the gapless and the gapped phases, the size correction 
decays sufficiently rapidly, so that A(g, L —oo) = ck. 
However, at the critical point, the estimator A(g, L)L/2'k 
instead converges to c -b B{q —>• 0)/27r. 


2. Paramagnetic phase 

For completeness we next extract the velocity of the 
propagating triplet excitations, or triplons, also in the 
quantum-disordered phase of the bilayer. We choose two 
points well away from g = at g = 3 and g = 4. The 
lowest triplet gap at (tt, tt, tt) now converges to a finite (g 
dependent) value as L —oo, because the paramagnetic 
phase is gapped. In the vicinity of (tt, tt, tt), the spectrum 
is expected to behave generically as 

A(q) = ^A2 + cV^Ao + ^, (56) 

where Aq denotes the triplet gap at wavevector (tt, 7r,7r), 
c is the velocity of the gapped triplons, and g is again 
measured relative to (tt, 7 r, 7 r). We use both the ap¬ 
proaches (ways of taking g —>• 0 and L —oo) discussed 
above to estimate the triplon velocity c and they give 
consistent results in the gapped phase; see Fig [50] for 
the analysis at g = 3. We obtain c = 1.973(4) and 
c = 2.159(6) for g = 3 and g = 4, respectively 

The simple estimator (1051) should give the correct ve¬ 
locity in the gapped phase, as discussed above, and we 
further check for consistency of the approach by also an¬ 
alyzing the gap at g = 47 r/L, in addition to the small¬ 
est momentum q = 2 tt/L (both extrapolations shown in 
Fig HOI upper panel). The same velocity is also obtained 
using these extrapolated gap values representing the limit 
L —>• oo before g —>• 0, as shown in Fig [501 bottom panel. 

For completeness, we also show the velocity estimate 
c(L) obtained from the winding numbers with the cubic 
criterion in Fig. |5TJ c(L) converges (exponentially fast 
for large L) to a finite value as L —>■ oo. However, this 
estimate gives an incorrect (higher) velocity in the para¬ 
magnetic phase. This is expected as the triplons are not 
linearly dispersing excitations [see Eq. (1551) ]. The error 
in the estimate increases with the distance of g from the 
critical point, which is because c is a regular function of 
(g — gc) from both sides and the winding number estima¬ 
tor does work at the critical point. 
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FIG. 20: (Color online) Velocity extraction for the bilayer 
at g = 3 (paramagnetic phase). The top panel shows the 
size convergence of velocity estimates based on two momenta 
close to the antiferromagnetic momentum {q being the de¬ 
viation from this momentum). The middle panel shows the 
system-size dependence of the gaps at several values of q. The 
bottom panel shows the infinite-size extrapolated velocity de¬ 
fined versus q, obtained from the fits in the middle panel and 
graphed versus q^. The 5 —^ 0 extrapolation by a polynomial 
(red curve) agrees with the different order of the limits taken 
in the top panel. The inset shows the dispersion relation. 



1/L 


FIG. 21: (Golor online) The behavior of the inverse tempera¬ 
ture /3* at which the cubic criterion is satisfied for the bilayer 
at two couplings in the paramagnetic regime. The method 
does not produce the correct velocity when g > Qc- The val¬ 
ues plus-minus one error bars obtained using the gap method 
is shown with the double horizontal lines. 

VII. QUANTUM-CRITICAL SCALING AT T > 0 

With c of the critical bilayer determined to high preci¬ 
sion, we have an opportunity to test in detail quantum- 
critical scaling predictions based on large-TV calculations 
for the 0{N) model, which for V = 3 should describe 
the universal critical behavior . The T = 0 quantum- 
critical point influences the behavior of the system in a 
wide “fan” in the (g, T) plane, with cross-overs to differ¬ 
ent low-T behaviors away from set by the spin stiffness 
(for g < gc) and the spin gap (for g > g^)- Here we will 
test results for the temperature dependence of the uni¬ 
form magnetic susceptibility x(T) (per unit cell), 

). ( 57 ) 

and the specific heat, 

C=dB(T), (58) 

where the internal energy E = {H) /per unit cell is 
computed with the bilayer Hamiltonian (l50ll . In the SSE 
method each configuration has a fixed z magnetization 
(since this is a conserved quantity) and x is evaluated di¬ 
rectly according to Eq. dSil). C can be computed using an 
exact estimator based on the fluctuations of the number 
of operators in the sampled operator strings. In practice, 
however, this estimator is very noisy at low temperatures 
and we will therefore instead analyze the energy, which is 
simply related to the average number of SSE operators. 

A. Results of large-N calculations 

The large-TV approximations give the following leading- 
order low-T forms of the susceptibility and the specific 
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heat at the critical coupling 


X'{T) = At, (59) 

TTC^ 

where the constant A = 1.0760, and 

C'iT) = (60) 

57rc^ 

where C(3) = 1.20206. Since we in practice compute and 
analyze the internal energy, we will compare with the 
T-integrated version of O; 


E'{T) =Eo + (61) 

We use primes on the symbols above to indicate that 
these are not expected to be exact forms. Several tests 
have been reported in the literature for different variants 
of dimerized Heisenberg models and the leading power 
laws above have been confirmed. However, quantitatively 
the degree of agreement has not been that well estab¬ 
lished, because c has also typically been extracted from 
quantities relying on the large-7V expansions, instead of 
using an independently determined unbiased value. 

One can also consider the Wilson ratio. 


W=^, (62) 

which has the advantage that its approximant based on 
Eqs. (15^ and (ISOl) does not involve the velocity. From 
the above expressions for % and C one obtains 


W' = 


5A 

12C(3) 


0.1243. 


(63) 


A value differing from this prediction by only about 2% 
was recently reported in Ref. [13 based on large-scale stud¬ 
ies of a single-plane model with columnar dimerization; 
W = 0.1262(6). Here our main aim is to obtain pre¬ 
cise estimates for the prefactors a and b of the leading 
low-T forms x = q’T, C = bT^ and compare these with 
the predicted values in with the value of the velocity 
c = 1.9001(2), as determined in the previous section. 


B. Susceptibility 

We need the susceptibility per unit cell in the thermo¬ 
dynamic limit. The finite-size effects are illustrated in 
Fig. for a coupling ratio close to the quantum-critical 
value. To make the graph clearer we have graphed x/T, 
which also is appropriate considering that we are here 
aiming to extract the size of the T-linear term in y. Since 
the magnetization is conserved and there is a finite size 
singlet-triplet gap A ~ 1/L, the susceptibility for a given 
finite L vanishes exponentially at a temperature scale 
T ~ 1/L. Interestingly x/T exhibits a prominent peak 
preceding the eventual drop to zero. The results indicate 



FIG. 22: (Color online) Susceptibility per unit cell divided 
by the temperature of the bilayer at ^ = 2.5220. Result for 
several different system sizes are shown to illustrate the finite 
size effects. For any finite L the susceptibility vanishes when 
T —>■ 0, at a temperature scale set by the finite size gap A, 
which at criticality scales as A ~ l/L. The lines connecting 
points are guides to the eye. 



FIG. 23: (Golor online) The susceptibility per spin of the 
bilayer Heisenberg model at coupling ratio s = 2.5221. The 
points are SSE results (with error bars much smaller than the 
symbols) and the curves are second-order polynomial fits to 
the data for T/Ji < 0.1. 


that we can compute the susceptibility without notice¬ 
able (within statistical errors) finite-size effects down to 
Tj J\ = 1/32 by using L = 512 lattices at the lowest tem¬ 
peratures (while smaller systems can be used at higher 
temperatures). 

Below we show results for a range of temperatures rep¬ 
resenting the thermodynamic limit. In addition to calcu¬ 
lations at g = 2.5220, we have also considered couplings 
one standard deviation of the gc estimate away from this 
point, i.e., g = 2.5219 and 2.5221. Based on these cal¬ 
culations we observe that the critical point should be 
very close to 2.5221. In Fig. [53] we analyze the T de¬ 
pendence of the susceptibility. At g = 2.5221 a second- 
order polynomial fit works well with data for Tj J\ < 0.1 
(x^/dof Ri 1.2 with 16 data points) and gives zero inter¬ 
cept within the statistical error of the fit; 0.000004(16). 

















21 



FIG. 24: (Color online) Twice the susceptibility divided by 
the temperature of the bilayer Heisenberg model at two cou¬ 
pling ratios close to the quantum-critical point. The solid 
curve is a third-order polynomial fit to the g = 2.5221 (which 
we estimate is closer to g^) data. 

The slope a in the expected leading-order form x = oT 
is a = 0.0922(3). If the intercept is assumed to be 0 and 
the form x = o-T bT^ is used, the resulting slope is 
0.09224(9). A second-order fit to data at at 5 = 2.5220 
and the same range of temperatures (not shown in the 
figure) gives a positive intercept about two error bars 
away from 0; 0.000032(14). The slope is a = 0.0916(6), 
not much different from the result for g = 2.5221. Al¬ 
though the differences between these data sets are small, 
they, along with results for g = 2.5219, suggest that the 
critical point is closer to g = 2.5221 than to i? = 2.5220, 
which is still in agreement within error bars with the 
result g— = 2.5220(1) stated in Ref. [13. Based on the 
present susceptibility results we estimate g = 2.52210(5). 

Fig. [M] shows the susceptibility per unit cell with the 
temperature divided out a.t g = 2.5220 and 2.5221, along 
with a fit to the data at the latter coupling. A third-order 
polynomial describe all the data well for all temperatures 
up to about T/Ji = 0.5. The intercept is completely con¬ 
sistent with the results for the slope obtained in Fig. [23l 
a = 0.09220(5). Given that this fit includes the largest 
amount of data we take the result as our final estimate of 
the susceptibility prefactor. We can now compare it with 
the predicted value from Eq. (I59L which with the value 
of c extracted in Sec. [Vl] is a' = 0.09487(2), i.e., 2.9% 
higher than our estimate. This close agreement is quite 
remarkable, considering that a' is based on a leading 1/N 
calculation evaluated at IV = 3. 


C. Specific heat 

The internal energy is graphed in Fig. 1251 The form 
E = Eg + bT^ describes well the data for T <2 and 
a fit gives E = -2.253040(1) and b = 0.2462(6). The 
predicted factor from Eq. (ICTl) is b' = 0.25435(5), which 
is 3.3% higher than our estimate. 



FIG. 25: (Color online) The internal energy of the bilayer 
Heisenberg model at coupling ratio g = 5.5221 (within the 
error bars of the estimated gEj, along with a fit of the form 
E — Eo + aT^ to the T < 0.2 data. The inset shows the data 
after subtraction of the ground-state energy and dividing by 
T®. The horizontal line shows the prefactor of the cubic cor¬ 
rection (with the line thickness corresponding approximately 
to ± one error bar) from the fit in the main graph. 

D. Wilson ratio 

The Wilson ratio with the exact prefactors a and b 
in the asymptotic forms x = ttT and E = Eq + bT^ 
{C = 3bT^) is FF = 3a/b. With the values of a and b 
determined above we obtain W = 0.1248(3). This value 
agrees reasonably well with the less accurate (with larger 
error bar) value obtained in Ref. being smaller by 
about two error bars. It is in remarkably good agreement 
with the 1/N value in Eq. (1551) . deviating by only one 
error bar. In other word, at the 95% confidence level 
(about two error bars), the 1/N estimate agrees with the 
true value to within about 0.5% or better. 


VIII. SUMMARY AND DISCUSSION 

We have here presented several methods to extract the 
velocity of excitations in quantum spin models simulated 
using QMC techniques. The methods are fundamentally 
not new, but we presented several technical solutions to 
improve their practical utility, and carried out extensive 
tests with the goal of testing their precision in practice. 

The method of computing the lowest excitation en¬ 
ergy for given momentum using the long-time behavior 
of imaginary-time correlations directly probes the dis¬ 
persion relation and is, thus, directly connected to the 
definition of the velocity. This is also the most compli¬ 
cated method in practice, as it requires significant ef¬ 
forts to obtain high-precision results for appropriate cor¬ 
relation functions and to analyze them, trying to avoid 
contamination by excited states in the gap estimation. 
We here employed two different methods to compute 
the correlations, using ground-state projector and finite- 
temperature QMC simulations. We also presented two 
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different ways to extract the gap using fitting procedures 
sensitive to the slowest decaying component of the cor¬ 
relation functions. The good agreement of the velocity 
obtained in this way with the other approaches in the 
limit of large system sizes show that the gap estimators 
are unbiased (and we also demonstrated this directly for 
small system sizes where exact results can be obtained). 
We expect these gap-extraction procedures to be useful 
also for many other systems and studies. 


The method of using winding number fluctuations to 
define a cubic space-time geometry fundamentally relies 
on the low-energy physics of systems with dynamic expo¬ 
nent z = 1. The method requires calculations for several 
temperatures in the neighborhood of the target tempera¬ 
ture T* at which the temporal and spatial winding num¬ 
ber fluctuations are equal, or, in principle, one could use 
a re-weighting technique and only do a single simulation 
very close to the estimated crossing point, though we ex¬ 
pect that it should be statistically advantageous to do 
several independent runs as more independent data are 
collected and the statistical analysis (curve fitting and er¬ 
ror estimate) is much simpler. Though doing several sim¬ 
ulations may seem like a drawback, in practice the pro¬ 
cedure for extracting the crossing temperature are very 
straight-forward and based on our work presented here 
this is the preferred method to extract the velocity. As 
noted in a previous application of this method,— the con¬ 
vergence with the system size is very rapid, but here we 
still found it necessary to use a final extrapolation using 
a power-law correction (with integer or non-trivial lead¬ 
ing power in a long-range ordered and critical system, 
respectively), to avoid any remaining size effects. 


We also presented a modification of the standard hy¬ 
drodynamic relationship between the velocity, spin stiff¬ 
ness, and susceptibility, computing the latter at non-zero 
momentum q = 2 tt/L in order to obtain a useful finite- 
size estimate of the velocity (with the standard relation¬ 
ship at q = 0 diverging when T —>■ 0 in a finite system 
in QMC simulations). This method also exhibits good 
convergence properties, though in practice it is still not 
as powerful as the winding-number method. 
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Appendix A: Nature of the systematical errors in 
the (n, dint) gap estimator 


We will show analytically the systematic error of the 
gap estimator ill) and the limits (1381) in more details, 
which was shown in Sec. Ill Bl Let us recall the dynamical 
correlation function in the exponential-form projection; 






where Z' = {-iptle IV't) = E, Cg|g), H\q) = Eg\q), 

and be^p^q = CpCq{p\A^C){^\A\q). The Fourier trans¬ 
formed correlation function at frequency ujk = 27rfc//3int 
ik G Z), is expressed as 


^ r^ivLt 

C{uJk)= / drC'(T,ro)e*' 

Jo 


(Al) 


1 

X] ^ X] di^pPint} (wfc = 0), 

Ee=Ep 


Finally, we extracted a high-precision estimate of the 
velocity in a critical bilayer model and used it in combi¬ 
nation with QMC calculated thermodynamic properties 
to test field-theory predictions based on 1 /N expansions 
of the 2-1-1 dimensional 0(3) model. The independently 
extracted velocity allowed us to precisely determine the 
deviations of the predicted overall factors in the temper¬ 
ature dependence of the susceptibility and the specific 
heat; the 1/A results (with N = 3) deviate from the true 
values by about 3%. The Wilson ratio, where the veloc¬ 
ity dependence cancels out, agrees with the 1/A value to 
within our 0.3% error bar (and we cannot establish the 
level of agreement beyond the level of the error bar). 


where Ant < 13 - tq, ge^p = de^p{l - e di^p = 

Ai^p = Ee - Ep. Note that 
the imaginary part for ujk ^ 0 does not vanish as 
finite-temperature case because the correlation function 
C{t,tq) is not periodic. 

The gap estimator (1571) is rewritten as 


_ 1 + An (Ant) 

Al 1 + Fn(Ant) + T)n(Ant) 



(A2) 


(Ant t OO ), 
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where 


^n(Ant) 


K(Ant) 


DniPint) 



h{n,l3int,£,p) 

Anti li 0) 



h{n,PintA,p) 
h(^Tl^ Anti li 0) 




Ant Ai 


1 

/l(n, Anti liO) 


, 

h{n, Pint, i,p) = (-l)”a;J” TT 2 , 2 ’ 

k=i ^e,p+^k 

be = be, 0,0 = |coA|^|0)p, and Ae = Eg - Eo- For the 
summation of i?„ and Fn, the set S is such that Ee ^ Ep 
except {i,p) = (liO). Then 


lim lim 

n-^oo /3i„t,To-too 


A 


("iftnt) 


= Ai. 


(A3) 


Let us consider taking the limit n —>■ oo before Anti To 
oo. Using the product expansion form of the hyperbolic 
function, sinh( 7 r 2 ) = + z'^jk'^), we obtain the 

limiting form as 

A(n,ft„t) l + -Roo(Ant) 10 

Ai 1 + Goo (Ant) 



where 


-Roo(Ant) — 


E 

Rples 






Goo (Ant) — 



,-%i(A,,p-Ai) 


(A5) 


Note that /3 —tq —since /3 > ro + Ant, Ap,o > 0 
(p 0 ), and Ae,i > 0 {£ > 1 ). Similarly, i?oo(Antj 0 
(Ant —t oo). Consequently, 


lim lim = Ai. (A 8 ) 

ftnt-S-OO n->-oo 

The convergence of Ant —>■ oo limit is accelerated together 
with taking tq —>■ oo limit practically. Therefore we have 
analytically shown that the limits are interchangeable: 


lim lim 

n-).oo Pint,To-too 


A 


(Tl,/3int) 


lim lim A 

Sint,To-too n^OO 


(tl,/3int) 


= Ai. 


This appealing property makes the estimation robust in 
practice and allows for flexibility in how the limit is taken. 

As we have shown in the main text, we extrapolate the 
final gap estimation for n ^ oo first and for Ant oo 
later. It is because that data for large Ant will be noisy, 
and it is easier to take the limit for n —> oo first. Nonethe¬ 
less, we could invert the limits. Then let us consider the 
finite n correction for making the extrapolation more re¬ 
liable. Let < 5(^1 z) = nfe=i(l + 2 ^/fc^). We have already 
used the limit Q{n,z) sinh( 7 r 2 ;)/ 7 r 2 (n —>• oo). The 
finite n correction is expressed as 


sinh(7rz) 

Q[n, z) ^ -—— exp 


sinh(7rz) 


TTZ 


1 - 


: -t- 1 2(?7. -t- l)(n' -t- 2) 


nz 


:-kl 2{n+l){n + 2) 2{n+iy 


Here the asymptotic expansion of the Riemann zeta func¬ 
tion was used; 


E 


fe=i 


1 


1 


6 n+1 


1 

2(71 -|- 1)(71 -|- 2) 


(A9) 


— ^ 1 , 0 , and Goo (Ant) — fim^_^oo(Fn(Ant) T Llyj(Ant))- 
Let us next consider the limit Ant 00 . We can rewrite 
the correction term as 

Goo(Ant) = ^ 5^_^g-(/3-T„)A„o-%l(A,„-Ai) 

(Ap)#(1,0) 

(bp)#(l,0) 

0 (Ant 00 ), (A6) 


where 


^he,p,qe 



q 


(A7) 


As a result, the asymptotic systematic error of the gap 
estimator is rewritten for Ant, tq ^ 1 as 

^(n.Ant) "" + E (/) (AlO) 

t.>i k 1 / 

'i + ^+ (JL + J— 

n + \ 2(n-I-1) \n-|-l n-I-2 

where be > 0 {£> 1), Ae,i = Ae — Ai >0 and z^ = 

— Af) > 0 {£ > 1). In practice, we can use 

a fitting function /(n,Ant) = Ai -|- aexp (—6 Ant)(l + 
cjn d/v?) with positive real parameters a, b, c, and real 
parameter d to extrapolate the first gap, as demonstrated 
in Fig. [3 and Fig. [ 6 ] in the main text. (Note that the 
parameter d may be negative for small Ant-) 
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